###################################################################################################
# Hurricane Sandy Analysis ########################################################################
# Heersink, Jenkins, and Olson ####################################################################
###################################################################################################

############################################
# dataset without proprietary SHELDUS data #
############################################

  dat <- read.csv("sandy_dat_nosheldus.csv")

##########################
# hurricane sandy data ###
##########################

# this is a date-location call to 
# SHELDUS 18.1 for all US counties from 10-26 to 11-02, 2012

  sandy <- read.csv("./sandy_organization/sheldus_1026_112.csv")

  sandy <- sandy[,c("State.Name","County.Name","County.FIPS","CropDmg.ADJ.2012.","PropertyDmg.ADJ.2012.","CropDmgPerCapita.ADJ.2012.","PropertyDmgPerCapita.ADJ.2012.")]
  colnames(sandy) <- c("statename","county","fips","crop","property","crop_pc","property_pc")
  sandy$fips <- gsub("'","",fixed=T,sandy$fips)
  state_merger <- as.data.frame(cbind(c(state.abb,"DC"),toupper(c(state.name,"district of columbia")))); colnames(state_merger) <- c("state","statename")
  sandy <- merge(sandy,state_merger,by=c("statename"),all.x=T)

  sandy$county <- gsub(" ","",fixed=T,tolower(sandy$county))
  sandy$tot_pc <- sandy$crop_pc+sandy$property_pc
  sandy <- sandy[,c("state","fips","tot_pc")]

# median affected county

  affected <- sandy[!is.na(sandy$tot_pc) & sandy$tot_pc!=0,]
  
  median(affected$tot_pc)*10000
  
# number of affected counties
  
  nrow(sandy[!is.na(sandy$tot_pc) & sandy$tot_pc>0,])
  
# merge into main dataset
  
  dat <- merge(dat,sandy,by=c("fips","state"),all.x=T)
  dat[is.na(dat$tot_pc),"tot_pc"] <- 0
  dat$tot_use <- log(dat$tot_pc*10000+1)
  dat$binary <- ifelse(dat$tot_use>0,1,0)

###################################################################################################
# summary statistics ##############################################################################
################################################################################################### 
  
  stargazer(dat[,c("obama12","obama08","tot_use","copart","swing",
                   "income_change","unemployment_change","housing_change",
                   "damage_pc","counter")],
            summary=T,
            covariate.labels = c("Obama 2012 Two-Party Share",
                                 "Obama 2008 Two-Party Share",
                                 "ln(Sandy Damage per 10,000)",
                                 "Co-Partisan County",
                                 "Swing County",
                                 "Income Change, 1 Year",
                                 "Unemployment Change, 1 Year",
                                 "Home Price Change, 2 Years",
                                 "ln(Disaster Damage per 10,000), 2 Years",
                                 "Disaster Declarations"),
            summary.stat=c("n","mean","median","sd","min","max"),
            out="./results/sandy_sumstats.tex",
            title="Summary Statistics, Hurricane Sandy Analysis",
            label="sandy_sumstats",
            font.size = "footnotesize")

###################################################################################################
# analysis ########################################################################################
###################################################################################################

# in-text table

  col1_no <- felm(obama12~tot_use
                  + obama08
                  |0|0|fips,data=dat)
  
  col2_no <- felm(obama12~tot_use*(copart+swing) 
                  + obama08|0|0|fips,data=dat)
  
  col4_no <- felm(obama12~tot_use*(copart+swing)
                  +income_change+unemployment_change+housing_change+
                    damage_pc+counter
                  + obama08|0|0|fips,data=dat)
  
  col1_div <- felm(obama12~tot_use
                   + obama08
                   |Division|0|fips,data=dat)
  
  col2_div <- felm(obama12~tot_use*(copart+swing) 
                   + obama08|Division|0|fips,data=dat)
  
  col4_div <- felm(obama12~tot_use*(copart+swing)
                   +income_change+unemployment_change+housing_change+
                     damage_pc+counter
                   + obama08|Division|0|fips,data=dat)
  
  sandy_new_sg <- stargazer(col1_no,col1_div,col2_no,col2_div,col4_no,col4_div,
                            dep.var.labels = c("Obama Vote Share 2012"),
                            title="Table 1: Impact of Hurricane Sandy on Obama Vote Share",
                            label="sandy_new",
                            order=c(1,10,11,2,3,4,5,6,7,8,9,12),
                            covariate.labels = c("Damage",
                                                 "Co-Partisan X$ Damage",
                                                 "Swing X Damage",
                                                 "Copartisan",
                                                 "Swing",
                                                 "Income Change, 1 Year",
                                                 "Unemployment Change, 1 Year",
                                                 "Home Price Change, 2 Years",
                                                 "Disaster Damage, 2 Years",
                                                 "Declarations, 2 Years",
                                                 "Obama Vote Share, 2008",
                                                 "Constant"),
                            keep.stat = c("n","rsq"),no.space = T,
                            type="html",
                            add.lines = list(c("Census Division FE","","</sup>&check;</sup>","","</sup>&check;</sup>","","</sup>&check;</sup>")),
                            star.cutoffs = c(0.10,0.05),star.char = c("*","**"),
                            notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                            table.layout ="-ld-#-t-as-n",out="./results/sandy_new.html",
                            notes="Note: Entries are linear regression coefficients with robust standard errors in parentheses. <sup>*</sup>p<0.10, <sup>**</sup><0.05 (two-tailed test).")

# effect in median affected county

  median <- median(dat$tot_use[dat$tot_use>0])
  
  median_treat_copart <- sum(coef(col4_div)[c("tot_use","tot_use:copart")])*median
  median_treat_contra <- sum(coef(col4_div)[c("tot_use")])*median
  
  print(median_treat_copart-median_treat_contra)

# make figure

  me_contra <- summary(col4_div)$coefficients["tot_use",1:2]
  
  me_copart_coef <- summary(col4_div)$coefficients["tot_use",1]+summary(col4_div)$coefficients["tot_use:copart",1]
  me_copart_se <- sqrt(summary(col4_div)$coefficients["tot_use",2]^2+summary(col4_div)$coefficients["tot_use:copart",2]^2
                       +2*col4_div$clustervcv["tot_use","tot_use:copart"])
  me_copart <- c(me_copart_coef,me_copart_se)
  
  me_swing_coef <- summary(col4_div)$coefficients["tot_use",1]+summary(col4_div)$coefficients["tot_use:swing",1]
  me_swing_se <- sqrt(summary(col4_div)$coefficients["tot_use",2]^2+summary(col4_div)$coefficients["tot_use:swing",2]^2
                      +2*col4_div$clustervcv["tot_use","tot_use:swing"])
  me_swing <- c(me_swing_coef,me_swing_se)

  
  plot_dat <- as.data.frame(rbind(me_contra,me_copart,me_swing,summary(col4_div)$coefficients["tot_use:copart",1:2]))
  plot_dat$Estimand <- c("Marginal Effect, Contra-Partisan",
                         "Marginal Effect, Co-Partisan",
                         "Marginal Effect, Swing",
                         "Difference, Co-Partisan and Contra-Partisan")
  plot_dat$Estimand <- factor(plot_dat$Estimand,
                              levels=c("Marginal Effect, Contra-Partisan",
                                       "Marginal Effect, Swing",
                                       "Marginal Effect, Co-Partisan",
                                       "Difference, Co-Partisan and Contra-Partisan"))
  
  cairo_ps("./results/fig2.eps",
      width=10,height=5,fallback_resolution = 1200)
  print(
    ggplot(data=plot_dat,aes(x=Estimand,y=Estimate,
                             ymin=Estimate-2*`Cluster s.e.`,
                             ymax=Estimate+2*`Cluster s.e.`))+
      geom_point(size=4)+
      geom_errorbar(width=0.25,size=2)+
      geom_hline(yintercept = 0,colour="indianred",size=2,linetype=2)+
      theme_minimal()
  )
  dev.off()
  
  png("./results/fig2.png",
           width=10,height=5,res=1200,units="in")
  print(
    ggplot(data=plot_dat,aes(x=Estimand,y=Estimate,
                             ymin=Estimate-2*`Cluster s.e.`,
                             ymax=Estimate+2*`Cluster s.e.`))+
      geom_point(size=4)+
      geom_errorbar(width=0.25,size=2)+
      geom_hline(yintercept = 0,colour="indianred",size=2,linetype=2)+
      theme_minimal()
  )
  dev.off()


###################################################################################################
# robustness ######################################################################################
################################################################################################### 

# 5-part split in democratic vote

  dat$quantile <- NA
  dat[dat$obama08<35,"quantile"] <- 1
  dat[dat$obama08>=35 & dat$obama08<45,"quantile"] <- 2
  dat[dat$obama08>=45 & dat$obama08<=55,"quantile"] <- 3
  dat[dat$obama08>55 & dat$obama08<=65,"quantile"] <- 4
  dat[dat$obama08>65,"quantile"] <- 5
  
  col2_quantile_no <- felm(obama12~tot_use*factor(quantile)
                           +obama08|0|0|fips,data=dat)
  
  col4_quantile_no <- felm(obama12~tot_use*factor(quantile)
                           +income_change+unemployment_change+housing_change
                           +damage_pc+counter
                           +obama08
                           |0|0|fips,data=dat)
  
  col2_quantile_div <- felm(obama12~tot_use*factor(quantile)
                            +obama08|Division|0|fips,data=dat)
  
  col4_quantile_div <- felm(obama12~tot_use*factor(quantile)
                            +income_change+unemployment_change+housing_change
                            +damage_pc+counter
                            +obama08
                            |Division|0|fips,data=dat)
  
  sandy_quantile_sg <- stargazer(col2_quantile_no,col2_quantile_div,col4_quantile_no,col4_quantile_div,
                                 dep.var.labels = c("Obama Vote Share 2012"),
                                 title="Impact of Sandy on Obama Vote Share: 5-Bin Measure of Partisanship",
                                 label="sandy_quantile",
                                 order=c(1,12,13,14,15,2,3,4,5,6,7,8,9,10,11),
                                 covariate.labels = c("Damage",
                                                      "Damage $\\times$ 35-45\\%",
                                                      "Damage $\\times$ 45-55\\%",
                                                      "Damage $\\times$ 55-65\\%",
                                                      "Damage $\\times$ >65\\%",
                                                      "35-45\\% for Obama",
                                                      "45-55\\% for Obama",
                                                      "55-65\\% for Obama",
                                                      ">65\\% for Obama",
                                                      "Income Change, 1 Year",
                                                      "Unemployment Change, 1 Year",
                                                      "Home Price Change, 2 Years",
                                                      "Disaster Damage, 2 Years",
                                                      "Declarations, 2 Years",
                                                      "Obama 2008 Vote Share",
                                                      "Constant"),
                                 add.lines=list(c("Census Division FE","","\\checkmark","","\\checkmark")),
                                 keep.stat = c("n","rsq"),no.space = T,
                                 star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!htbp",
                                 notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                                 table.layout ="-ld-#-t-as-n",font.size="footnotesize",
                                 notes="\\parbox[t]{0.75\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       robust standard errors in parentheses. $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(sandy_quantile_sg, sep = '\n', file = "./results/sandy_quantile.tex")

# make figure

  me_1 <- summary(col4_quantile_div)$coefficients["tot_use",1:2]
  
  me_2_coef <- summary(col4_quantile_div)$coefficients["tot_use",1]+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)2",1]
  me_2_se <- sqrt(summary(col4_quantile_div)$coefficients["tot_use",2]^2+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)2",2]^2
                  +2*col4_quantile_div$clustervcv["tot_use","tot_use:factor(quantile)2"])
  me_2 <- c(me_2_coef,me_2_se)
  
  me_3_coef <- summary(col4_quantile_div)$coefficients["tot_use",1]+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)3",1]
  me_3_se <- sqrt(summary(col4_quantile_div)$coefficients["tot_use",2]^2+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)3",2]^2
                  +2*col4_quantile_div$clustervcv["tot_use","tot_use:factor(quantile)3"])
  me_3 <- c(me_3_coef,me_3_se)
  
  me_4_coef <- summary(col4_quantile_div)$coefficients["tot_use",1]+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)4",1]
  me_4_se <- sqrt(summary(col4_quantile_div)$coefficients["tot_use",2]^2+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)4",2]^2
                  +2*col4_quantile_div$clustervcv["tot_use","tot_use:factor(quantile)4"])
  me_4 <- c(me_4_coef,me_4_se)
  
  me_5_coef <- summary(col4_quantile_div)$coefficients["tot_use",1]+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)5",1]
  me_5_se <- sqrt(summary(col4_quantile_div)$coefficients["tot_use",2]^2+summary(col4_quantile_div)$coefficients["tot_use:factor(quantile)5",2]^2
                  +2*col4_quantile_div$clustervcv["tot_use","tot_use:factor(quantile)5"])
  me_5 <- c(me_5_coef,me_5_se)
  
  plot_dat <- as.data.frame(rbind(me_1,me_2,me_3,me_4,me_5))
  plot_dat$Estimand <- c("ME, <35% for Obama",
                         "ME, 35-45% for Obama",
                         "ME, 45-55% for Obama",
                         "ME, 55-65% for Obama",
                         "ME, >65% for Obama")
  plot_dat$Estimand <- factor(plot_dat$Estimand,
                              levels=c("ME, <35% for Obama",
                                       "ME, 35-45% for Obama",
                                       "ME, 45-55% for Obama",
                                       "ME, 55-65% for Obama",
                                       "ME, >65% for Obama"))
  
  pdf("./results/quantile_plot.pdf",width=10,height=5)
  print(
    ggplot(data=plot_dat,aes(x=Estimand,y=Estimate,
                             ymin=Estimate-2*`Cluster s.e.`,
                             ymax=Estimate+2*`Cluster s.e.`))+
      geom_point(size=4)+
      geom_errorbar(width=0.25,size=2)+
      geom_hline(yintercept = 0,colour="indianred",size=2,linetype=2)+
      theme_minimal()
  )
  dev.off()

# Continuous vote share measure  

  col2_cont_no <- felm(obama12~tot_use*obama08|0|0|fips,data=dat)
  
  col4_cont_no <- felm(obama12~tot_use*obama08
                       +income_change+unemployment_change+housing_change
                       +damage_pc+counter
                       |0|0|fips,data=dat)
  
  col2_cont_div <- felm(obama12~tot_use*obama08|Division|0|fips,data=dat)
  
  col4_cont_div <- felm(obama12~tot_use*obama08
                        +income_change+unemployment_change+housing_change
                        +damage_pc+counter
                        |Division|0|fips,data=dat)
  
  sandy_cont_sg <- stargazer(col2_cont_no,col2_cont_div,col4_cont_no,col4_cont_div,
                             dep.var.labels = c("Obama Vote Share 2012"),
                             title="Impact of Hurricane Sandy on Obama Vote Share: Continuous Measure of Partisanship",
                             label="sandy_cont",
                             order=c(1,2,8,3,4,5,6,7,9),
                             covariate.labels = c("Damage",
                                                  "Obama Vote Share, 2008",
                                                  "Damage $\\times$ Obama Vote Share, 2008",
                                                  "Income Change, 1 Year",
                                                  "Unemployment Change, 1 Year",
                                                  "Home Price Change, 2 Years",
                                                  "Disaster Damage, 2 Years",
                                                  "Declarations, 2 Years",
                                                  "Constant"),
                             add.lines = list(c("Census Division FE","","\\checkmark","","\\checkmark")),
                             keep.stat = c("n","rsq"),no.space = T,
                             star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                             notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                             table.layout ="-ld-#-t-as-n",
                             notes="\\parbox[t]{0.85\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       robust standard errors in parentheses. $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(sandy_cont_sg, sep = '\n', file = "./results/sandy_cont.tex")

# present this as a marginal effects plot, too

  cont_me <- inter.kernel(data=dat,Y="obama12",D="tot_use",X="obama08",
                          Z=c("income_change","unemployment_change","housing_change",
                              "damage_pc","counter"),
                          treat.type="continuous",
                          full.moderate=TRUE,FE="Division",
                          na.rm=T,CI=T,main="",ylab="Marginal Effect of Sandy on Obama 2012 Vote Share",
                          xlab="Obama 2008 Vote Share",theme.bw = T)
  
  pdf("./results/cont_plot.pdf",width=10,height=6.5)
  print(cont_me$graph)
  dev.off()

# binary measure of sandy damage

  col1_bin_no <- felm(obama12~binary
                      + obama08
                      |0|0|fips,data=dat)
  
  col2_bin_no  <- felm(obama12~binary*(copart+swing) 
                       + obama08|0|0|fips,data=dat)
  
  col4_bin_no  <- felm(obama12~binary*(copart+swing) 
                       +income_change+unemployment_change+housing_change
                       +damage_pc+counter
                       + obama08|0|0|fips,data=dat)
  
  col1_bin_div <- felm(obama12~binary
                       + obama08
                       |Division|0|fips,data=dat)
  
  col2_bin_div  <- felm(obama12~binary*(copart+swing) 
                        + obama08|Division|0|fips,data=dat)
  
  col4_bin_div  <- felm(obama12~binary*(copart+swing) 
                        +income_change+unemployment_change+housing_change
                        +damage_pc+counter
                        + obama08|Division|0|fips,data=dat)
  
  sandy_bin_new_sg <- stargazer(col1_bin_no ,col1_bin_div ,col2_bin_no ,
                                col2_bin_div ,col4_bin_no ,col4_bin_div ,
                                dep.var.labels = c("Obama Vote Share 2012"),
                                title="Impact of Hurricane Sandy on Obama Vote Share: Binary Measure of Sandy Damage",
                                label="sandy_bin",
                                order=c(1,10,11,2,3,4,5,6,7,8,9,12),
                                covariate.labels = c("Sandy Dummy",
                                                     "Co-Partisan $\\times$ Sandy",
                                                     "Swing $\\times$ Sandy",
                                                     "Co-Partisan",
                                                     "Swing",
                                                     "Income Change, 1 Year",
                                                     "Unemployment Change, 1 Year",
                                                     "Home Price Change, 2 Years",
                                                     "Disaster Damage, 2 Years",
                                                     "Declarations, 2 Years",
                                                     "Obama Vote Share, 2008",
                                                     "Constant"),
                                add.lines=list(c("Census Division FE","","\\checkmark","","\\checkmark","","\\checkmark")),
                                keep.stat = c("n","rsq"),no.space = T,font.size="footnotesize",
                                star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                                notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                                table.layout ="-ld-#-t-as-n",
                                notes="\\parbox[t]{0.95\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       robust standard errors in parentheses. $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(sandy_bin_new_sg, sep = '\n', file = "./results/sandy_bin_new.tex")

# alternative control variable specifications

  col4_salary_no <- felm(obama12~tot_use*(copart+swing)
                         +income_change_alt+unemployment_change+housing_change+
                           damage_pc+counter
                         + obama08|0|0|fips,data=dat)
  
  col4_inc4_no <- felm(obama12~tot_use*(copart+swing)
                       +income_change08+unemployment_change+housing_change+
                         damage_pc+counter
                       + obama08|0|0|fips,data=dat)
  
  col4_housing4_no <- felm(obama12~tot_use*(copart+swing)
                           +income_change+unemployment_change+housing_change4+
                             damage_pc+counter
                           + obama08|0|0|fips,data=dat)
  
  col4_noem_no <- felm(obama12~tot_use*(copart+swing)
                       +income_change+unemployment_change+housing_change+
                         damage_pc+counter_noem
                       + obama08|0|0|fips,data=dat)
  
  col4_salary_div <- felm(obama12~tot_use*(copart+swing)
                          +income_change_alt+unemployment_change+housing_change+
                            damage_pc+counter
                          + obama08|Division|0|fips,data=dat)
  
  col4_inc4_div <- felm(obama12~tot_use*(copart+swing)
                        +income_change08+unemployment_change+housing_change+
                          damage_pc+counter
                        + obama08|Division|0|fips,data=dat)
  
  col4_housing4_div <- felm(obama12~tot_use*(copart+swing)
                            +income_change+unemployment_change+housing_change4+
                              damage_pc+counter
                            + obama08|Division|0|fips,data=dat)
  
  col4_noem_div <- felm(obama12~tot_use*(copart+swing)
                        +income_change+unemployment_change+housing_change+
                          damage_pc+counter_noem
                        + obama08|Division|0|fips,data=dat)
  
  sandy_alt_new_sg <- stargazer(col4_salary_no,col4_salary_div,col4_inc4_no,col4_inc4_div,
                                col4_housing4_no,col4_housing4_div,col4_noem_no,col4_noem_div,
                                dep.var.labels = c("Obama Vote Share 2012"),
                                title="Impact of Hurricane Sandy on Obama Vote Share: Alternative Control Variables",
                                label="sandy_alt_new",
                                order=c(1,14,15,2,3,4,6,5,7,8,9,10,11,12,13,16),
                                covariate.labels = c("Damage",
                                                     "Co-Partisan $\\times$ Damage",
                                                     "Swing $\\times$ Damage",
                                                     "Co-Partisan",
                                                     "Swing",
                                                     "$\\Delta$ Salary, 1 Year",
                                                     "$\\Delta$ Income, 1 Year",
                                                     "$\\Delta$ Income, 4 Years",
                                                     "$\\Delta$ Unemp., 1 Year",
                                                     "$\\Delta$ Home Price, 2 Years",
                                                     "$\\Delta$ Home Price, 4 Years",
                                                     "Disaster Damage, 2 Years",
                                                     "Declarations, 2 Years",
                                                     "Declarations, No Emerg.",
                                                     "Obama Vote Share, 2008",
                                                     "Constant"),
                                add.lines=list(c("Census Division FE","","\\checkmark","","\\checkmark","","\\checkmark","","\\checkmark")),
                                keep.stat = c("n","rsq"),no.space = T,font.size="footnotesize",
                                star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                                notes.append = FALSE,notes.label = "",column.sep.width="-5pt",
                                table.layout ="-ld-#-t-as-n",
                                notes="\\parbox[t]{\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       robust standard errors in parentheses. $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(sandy_alt_new_sg, sep = '\n', file = "./results/sandy_alt_new.tex")

# controlling for additional lags and alternative fixed effects

  col4_3lag_no  <- felm(obama12~tot_use*(copart+swing) 
                        +income_change+unemployment_change+housing_change
                        +damage_pc+counter
                        + obama08 + kerry04 + gore00|0|0|fips,data=dat)
  
  col4_3lag_div  <- felm(obama12~tot_use*(copart+swing) 
                         +income_change+unemployment_change+housing_change
                         +damage_pc+counter
                         + obama08 + kerry04 + gore00|Division|0|fips,data=dat)
  
  col4_state <- felm(obama12~tot_use*(copart+swing)
                    +income_change+unemployment_change+housing_change
                    +damage_pc+counter
                     + obama08
                     |state|0|fips,data=dat)
  
  col4_region <- felm(obama12~tot_use*(copart+swing)
                      +income_change+unemployment_change+housing_change
                      +damage_pc+counter
                      + obama08
                      |Region|0|fips,data=dat)
  
  sandy_lag_new_sg <- stargazer(col4_3lag_no,col4_3lag_div,col4_state,col4_region,
                                dep.var.labels = c("Obama Vote Share 2012"),
                                title="Impact of Hurricane Sandy on Obama Vote Share: Alternative Lags and Fixed Effects",
                                label="sandy_lag_new",
                                order=c(1,12,13,2,3,4,5,6,7,8,9,10,11,14),
                                covariate.labels = c("Damage",
                                                     "Co-Partisan $\\times$ Damage",
                                                     "Swing $\\times$ Damage",
                                                     "Co-Partisan",
                                                     "Swing",
                                                     "Income Change, 1 Year",
                                                     "Unemployment Change, 1 Year",
                                                     "Home Price Change, 2 Years",
                                                     "Disaster Damage, 2 Years",
                                                     "Declarations, 2 Years",
                                                     "Obama Vote Share, 2008",
                                                     "Kerry Vote Share, 2004",
                                                     "Gore Vote Share, 2000",
                                                     "Constant"),
                                add.lines = list(c("Fixed Effects","","Division","State","Census Region")),
                                keep.stat = c("n","rsq"),no.space = T,
                                star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                                notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                                table.layout ="-ld-#-t-as-n",font.size="footnotesize",
                                notes="\\parbox[t]{0.75\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       robust standard errors in parentheses. $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(sandy_lag_new_sg, sep = '\n', file = "./results/sandy_lag_new.tex")  

# fully moderated models

  col4_mod_no <- felm(obama12~(copart+swing)*(tot_use
                                              +income_change+unemployment_change+housing_change+
                                                damage_pc+counter
                                              + obama08)|0|0|fips,data=dat)
  
  col4_mod_div <-  felm(obama12~(copart+swing)*(tot_use
                                                +income_change+unemployment_change+housing_change+
                                                  damage_pc+counter
                                                + obama08+factor(Division))-factor(Division)|Division|0|fips,data=dat)
  
  sandy_mod_sg <- stargazer(col4_mod_no,col4_mod_div,
                            dep.var.labels = c("Obama Vote Share 2012"),
                            title="Impact of Hurricane Sandy on Obama Vote Share: Fully Moderated Models",
                            label="sandy_mod",
                            keep=c("tot_use"),
                            covariate.labels = c("Damage","Co-Partisan $\\times$ Damage","Swing $\\times$ Damage"),
                            add.lines = list(c("Census Division FE","","\\checkmark")),
                            keep.stat = c("n","rsq"),no.space = T,single.row = T,
                            star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                            notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                            table.layout ="-ld-#-t-as-n",
                            notes="\\parbox[t]{0.65\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       robust standard errors in parentheses. All models include all covariates from Model 4 in Table 1 and interactions between 
                       `Co-Partisan' and `Swing' and these covariates. Model 2 additionally includes interactions between Census Division FE and 
                       `Co-Partisan' and `Swing.' $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(sandy_mod_sg, sep = '\n', file = "./results/sandy_mod.tex")

# incorporating ideology  
  
  print(cor(dat$obama08,dat$mrp_ideology_mean,use="complete.obs"))
  

